��#  
 #  
 # R   C o d e   f o r   e s t i m a t i n g   Z I O P   a n d   Z I O P C   m o d e l s  
 #  
 # B E B   4 / 2 7 / 2 0 1 1  
 #  
 #  
 #  
  
 # c l e a r   m e m o r y    
 r m (   l i s t = l s ( )   )    
  
 # l o a d   t w o   n e c e s s a r y   p a c k a g e s  
 l i b r a r y ( m v t n o r m )  
 l i b r a r y ( M A S S )  
  
 # s e t   t h e   s e e d  
 s e t . s e e d ( 2 )  
  
 # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #  
 # # # # # # # # # # # C r e a t e   a   d a t a s e t   f o r   u s e   i n   t h e   m o d e l s   b e l o w # # # # # # # # # # #  
 # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #  
  
 # s e t   t h e   n u m b e r   o f   o b s e r v a t i o n s  
     n < - 1 0 0 0  
  
 # d e f i n e   t r u e   g a m m a s   a n d   b e t a s  
     G < - r b i n d ( 1 ,   - 0 . 2 5 ,   - 1 )        
     G < - a s . m a t r i x ( G )          
     c u t p o i n t s < - r b i n d ( - 1 0 0 ,   0 ,   4 . 5 ,   5 . 5 ,   1 0 0 )  
     B < - r b i n d ( 1 , 1 )  
     B < - a s . m a t r i x ( B )  
  
 # c r e a t e   s o m e   c o v a r i a t e s  
     x 0 < - m a t r i x ( r e p ( 1 , n ) )                                                                                                                        
     x 1 < - m a t r i x ( l o g ( r u n i f ( n ,   m i n = 0 ,   m a x = 1 0 0 ) ) , n , 1 )  
     x 2 s t a r < - m a t r i x ( r u n i f ( n ,   m i n = 0 ,   m a x = 1 ) ,   n ,   1 )              
     x 2 < - i f e l s e ( x 2 s t a r > . 2 5 , 1 , 0 )        
     Z < - c b i n d ( x 0 , x 1 , x 2 )  
     Z < - a s . m a t r i x ( Z )      
     X < - c b i n d ( x 1 , x 2 )  
     X < - a s . m a t r i x ( X )  
  
 # c r e a t e   a   z e r o   i n f l a t e d   d e p e n d e n t   v a r i a b l e   t h a t  
 # r a n g e s   f r o m   0   t o   3 ,   w i t h   n o   c o r r e l a t e d   d i s t u r b a n c e s  
     p h i < - p n o r m ( Z % * % G )    
     v c v < - m a t r i x ( c ( 1 , 0 , 0 , 1 ) ,   n r o w = 2 ,   n c o l = 2 )  
     e < - m v r n o r m ( n ,   m u = c ( 0 , 0 ) ,   S i g m a = v c v )  
     y s t a r < - X % * % B + e [ , 1 ]  
     y < - m a t r i x ( 0 , n , 1 )  
     f o r ( i   i n   2 : 4 ) { y [ y s t a r > c u t p o i n t s [ i ] ] < - y [ y s t a r > c u t p o i n t s [ i ] ] + 1 }  
     y z e r o < - m a t r i x ( 0 , n , 1 )  
     e r r o r < - - 1 * e [ , 2 ]  
     f l a g < - e r r o r < q n o r m ( p h i )  
     y z e r o [ f l a g ] < - e r r o r [ f l a g ]  
     f l a g < - y z e r o = = 0  
     y [ f l a g ] < - y z e r o [ f l a g ]  
  
 # s a v e   c o m p l e t e d   d a t a s e t  
     d a t a s e t < - c b i n d ( y , x 1 , x 2 )    
     c o l n a m e s ( d a t a s e t )   < -   c ( " y " ,   " x 1 " ,   " x 2 " )  
     d a t a s e t < - a s . d a t a . f r a m e ( d a t a s e t )  
  
  
 
 # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #  
 # # # # # # # # # # # # # # # # # # # # # # # # # # # # Z I O P   M o d e l # # # # # # # # # # # # # # # # # # # # # # # # # # # #  
 # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #  
  
  
 # T h i s   p r o g r a m   e s t i m a t e s   t h e   Z I O P   l o g l i k e l i h o o d   f u n c t i o n  
  
 Z I O P < -   f u n c t i o n ( e s t , Y , X , Z , d a t a )   { 	 	 	 	 	              
 	 n = n r o w ( d a t a ) 	 	 	 	 	 	 	             	 	 	 	 	      
 	 l l i k   < -   m a t r i x ( 0 ,   n r o w = n ,   n c o l   =   1 )  
 	 y . c a t < - n l e v e l s ( a s . f a c t o r ( Y ) )  
 	 y 0 < - s o r t ( u n i q u e ( Y ) )  
 	 V < - m a t r i x ( N A , n r o w ( d a t a ) , y . c a t )  
 	 f o r ( k   i n   1 : y . c a t ) {  
 	 	 V [ , k ] < - Y = = y 0 [ k ]  
 	 	 }  
 	 t a u < - r e p ( 1 , m a x ( Y ) )  
 	 f o r   ( i   i n   1 : m a x ( Y ) ) {  
 	 	 t a u [ i ] < - i f e l s e ( i = = 1 , e s t [ i ] , t a u [ i - 1 ] + e x p ( e s t [ i ] ) )  
 	 	 }  
 	 g a m m a < - e s t [ ( y . c a t ) : ( y . c a t + n c o l ( Z ) - 1 ) ]  
 	 b e t a < - e s t [ ( y . c a t + n c o l ( Z ) ) : l e n g t h ( e s t ) ]  
 	 Z G < - Z % * % g a m m a  
 	 X B < - X % * % b e t a 
 	  
 	 c p r o b s < - m a t r i x ( n r o w = l e n g t h ( X B ) , n c o l = y . c a t )  
 	 p r o b s < - m a t r i x ( n r o w = n , n c o l = y . c a t ) 
 	            
 	 f o r ( j   i n   1 : ( y . c a t - 1 ) ) { 
 	 	  
 	 	 c p r o b s [ , j ] < - p n o r m ( t a u [ j ] - X B ) 
 	 	  
 	 	 } 
 	  
 	 p r o b s [ , y . c a t ] < - ( 1 - c p r o b s [ , ( y . c a t - 1 ) ] ) * p n o r m ( Z G )    
 	 p r o b s [ , 1 ] < - ( c p r o b s [ , 1 ] ) * p n o r m ( Z G ) + ( 1 - p n o r m ( Z G ) ) 	 
 	  
 	 f o r ( j   i n   2 : ( y . c a t - 1 ) ) { 	 	 	 
 	 	  
 	 	 p r o b s [ , j ] < - ( c p r o b s [ , j ] - c p r o b s [ , ( j - 1 ) ] ) * ( p n o r m ( Z G ) ) 
 	 	  
 	 	 } 
 	  
 	 l l i k < - - 1 * s u m (   l o g ( p r o b s [ V ] )   ) 
 	  
 	 r e t u r n ( l l i k ) 
 	  
 	 } 
  
  
  
 # s e t   s t a r t i n g   p a r a m e t e r s  
 e s t < - r b i n d ( . 1 , . 1 , . 1 , . 1 , . 1 , . 1 , . 1 , . 1 )  
  
 # s e t   d a t a ,   Y ,   X ,   a n d   Z  
 d a t a < - d a t a s e t  
 Y < - d a t a s e t $ y  
 X < - c b i n d ( d a t a s e t $ x 1 , d a t a s e t $ x 2 )  
 Z < - c b i n d ( 1 ,   d a t a s e t $ x 1 ,   d a t a s e t $ x 2 )  
  
 # o p t i m i z e  
 o u t p u t . Z I O P < - n l m ( f = Z I O P ,     p = e s t ,   Y = Y , X = X , Z = Z   , i t e r l i m = 5 0 0 ,   d a t a = d a t a s e t ,   h e s s i a n = T R U E )  
 p r i n t ( o u t p u t . Z I O P )  
 v c v < - s o l v e ( o u t p u t . Z I O P $ h e s s i a n )  
 v c v    
  
  
  
  
 
 
 
 # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #  
 # # # # # # # # # # # # # # # # # # # # # # # # # # # # Z I O P C   M o d e l # # # # # # # # # # # # # # # # # # # # # # # # # # # # #  
 # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
 
 #  
  
 # T h i s   p r o g r a m   e s t i m a t e s   t h e   Z I O P C   l o g l i k e l i h o o d   f u n c t i o n  
  
 Z I O P C < -   f u n c t i o n ( e s t , Y , X , Z , d a t a )   { 	 	 	 	 	              
 	 n = n r o w ( d a t a ) 	 	 	 	 	 	 	             	 	 	 	 	      
 	 l l i k   < -   m a t r i x ( 0 ,   n r o w = n ,   n c o l   =   1 )  
 	 y . c a t < - n l e v e l s ( a s . f a c t o r ( Y ) )  
 	 y 0 < - s o r t ( u n i q u e ( Y ) )  
 	 V < - m a t r i x ( N A , n r o w ( d a t a ) , y . c a t )  
 	 f o r ( k   i n   1 : y . c a t ) {  
 	 	 V [ , k ] < - Y = = y 0 [ k ]  
 	 	 }  
 	 t a u < - r e p ( 1 , m a x ( Y ) )  
 	 f o r   ( i   i n   1 : m a x ( Y ) ) {  
 	 	 t a u [ i ] < - i f e l s e ( i = = 1 , e s t [ i ] , t a u [ i - 1 ] + e x p ( e s t [ i ] ) )  
 	 	 }  
 
 	 g a m m a < - e s t [ ( ( y . c a t ) ) : ( y . c a t + n c o l ( Z ) - 1 ) ]  
 	 b e t a < - e s t [ ( y . c a t + n c o l ( Z ) ) : ( l e n g t h ( e s t ) - 1 ) ]  
 	 r h o < - e s t [ l e n g t h ( e s t ) ] 
 	  
 	 s i g m a < - m a t r i x ( c ( 1 , r h o , r h o , 1 ) ,   n r o w = 2 ,   n c o l = 2 )  
 	 n e g s i g m a < - m a t r i x ( c ( 1 , - r h o , - r h o , 1 ) ,   n r o w = 2 ,   n c o l = 2 ) 
 	 
 	 
 	  
 	 c p r o b s < - m a t r i x ( n r o w = n , n c o l = y . c a t )  
 	 p r o b s < - m a t r i x ( n r o w = n , n c o l = y . c a t ) 
                
 	 f o r ( j   i n   1 : ( y . c a t - 1 ) ) { 
 	 	  
 	 	 f o r ( i   i n   1 : n ) { 
 	 	  
 	 	 	 Z G < - Z [ i , ] % * % g a m m a  
 	 	 	 X B < - X [ i , ] % * % b e t a 
 	 	  
 	 	 	 c p r o b s [ i , j ] < - p m v n o r m ( l o w e r = c ( - I n f , - I n f ) ,   u p p e r = c ( Z G , t a u [ j ] - X B ) , m e a n = c ( 0 , 0 ) , s i g m a = n e g s i g m a ) 	 
 	 	  
 	 	 	 } 	 
 	  
 	 	 } 	 
 	  
 	 f o r ( i   i n   1 : n ) { 
 	  
 	 	 Z G < - Z [ i , ] % * % g a m m a  
 	 	 X B < - X [ i , ] % * % b e t a 
 	  
 	 	 p r o b s [ i , 1 ] < - p m v n o r m ( l o w e r = c ( - I n f , - I n f ) ,   u p p e r = c ( Z G , t a u [ 1 ] - X B ) , m e a n = c ( 0 , 0 ) , s i g m a = n e g s i g m a ) + ( 1 - p n o r m ( Z G ) ) 
 	  
 	 	 } 
 	  
 	 f o r ( i   i n   1 : n ) { 	 
 	 	  
 	 	 Z G < - Z [ i , ] % * % g a m m a  
 	 	 X B < - X [ i , ] % * % b e t a 
 	 	  
 	 	 p r o b s [ i , y . c a t ] < - ( p m v n o r m ( l o w e r = c ( - I n f , - I n f ) ,   u p p e r = c ( Z G , ( X B - t a u [ y . c a t - 1 ] ) ) , m e a n = c ( 0 , 0 ) , s i g m a = s i g m a ) )     
 	 	  
 	 	 } 
 	  
 	 f o r ( j   i n   2 : ( y . c a t - 1 ) ) { 	 	 	 
 	 	  
 	 	 p r o b s [ , j ] < - ( c p r o b s [ , j ] - c p r o b s [ , ( j - 1 ) ] ) 
 	 	  
 	 	 } 
 	  
 	 l l i k < - - 1 * s u m (   l o g ( p r o b s [ V ] )   ) 
 	  
 	 r e t u r n ( l l i k ) 
 	  
 	 } 
  
  
 # s e t   s t a r t i n g   p a r a m e t e r s  
 e s t < - r b i n d ( . 1 , . 1 , . 1 , . 1 , . 1 , . 1 , . 1 , . 1 , . 1 )  
  
 # s e t   d a t a ,   Y ,   X ,   a n d   Z  
 d a t a < - d a t a s e t  
 Y < - d a t a s e t $ y  
 X < - c b i n d ( d a t a s e t $ x 1 , d a t a s e t $ x 2 )  
 Z < - c b i n d ( 1 ,   d a t a s e t $ x 1 ,   d a t a s e t $ x 2 )  
  
 # o p t i m i z e  
 o u t p u t . Z I O P C < - n l m ( f = Z I O P C ,     p = e s t ,   Y = Y , X = X , Z = Z   , i t e r l i m = 5 0 0 ,   d a t a = d a t a s e t ,   h e s s i a n = T R U E )  
 p r i n t ( o u t p u t . Z I O P C )  
 v c v < - s o l v e ( o u t p u t . Z I O P C $ h e s s i a n )  
 v c v    
  
  
  
 
 
 # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #  
 # # # # # # # # # # # # # # # # # # # A l t e r n a t e   Z I O P C   M o d e l # # # # # # # # # # # # # # # # # # # # # # # # # # # #  
 # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # 
 
 #  
  
 # T h i s   p r o g r a m   e s t i m a t e s   a   s l i g h t l y   f a s t e r   v e r s i o n   o f   t h e   Z I O P C   m o d e l  
  
 # F o r   t h i s   v e r s i o n ,   o n e   m u s t   s e t   t h e   n u m b e r   o f   Y   c a t e g o r i e s   m a n u a l l y  
 # ( h e r e   t h e   c a t e g o r i e s   a r e   s e t   t o   0 - 3 ,   t o   c o r r e s p o n d   t o   t h e   f o u r   c a t e g o r i e s   o f   o u r   Y ) .  
  
 Z I O P C   < -   f u n c t i o n ( e s t , Y , X , Z , d a t a )   { 	 	 	 	 	              
 	 n = n r o w ( d a t a ) 	 	 	 	 	 	 	             	 	 	 	 	 	 	              
 	 l l i k   < -   m a t r i x ( 0 ,   n r o w = n ,   n c o l   =   1 )  
 	 y . c a t < - n l e v e l s ( a s . f a c t o r ( Y ) )  
 	 t a u < - r e p ( 1 , m a x ( Y ) )  
 	 f o r   ( i   i n   1 : m a x ( Y ) ) {  
 	 	 t a u [ i ] < - i f e l s e ( i = = 1 , e s t [ i ] , t a u [ i - 1 ] + e x p ( e s t [ i ] ) )  
 	 	 }  
 	 g a m m a < - e s t [ ( ( y . c a t ) ) : ( y . c a t + n c o l ( Z ) - 1 ) ]  
 	 b e t a < - e s t [ ( y . c a t + n c o l ( Z ) ) : ( l e n g t h ( e s t ) - 1 ) ]  
 	 r h o < - e s t [ l e n g t h ( e s t ) ]  
 	 f o r ( i   i n   1 : n ) {  
 	 	 Z G < - Z [ i , ] % * % g a m m a  
 	 	 X B < - X [ i , ] % * % b e t a  
 	 	 s i g m a < - m a t r i x ( c ( 1 , r h o , r h o , 1 ) ,   n r o w = 2 ,   n c o l = 2 )  
 	 	 n e g s i g m a < - m a t r i x ( c ( 1 , - r h o , - r h o , 1 ) ,   n r o w = 2 ,   n c o l = 2 )  
 	 	 i f ( i d e n t i c a l ( Y [ i ] , 0 ) ) { l l i k [ i ] < - l o g ( ( 1 - p n o r m ( Z G ) ) +   p m v n o r m ( l o w e r = c ( - I n f , - I n f ) ,   u p p e r = c ( Z G , ( t a u [ 1 ] - X B ) ) , m e a n = c ( 0 , 0 ) , s i g m a = n e g s i g m a ) ) }  
 	 	 e l s e   i f ( i d e n t i c a l ( Y [ i ] , 1 ) ) { l l i k [ i ] < - l o g ( p m v n o r m ( l o w e r = c ( - I n f , - I n f ) ,   u p p e r = c ( Z G , ( t a u [ 2 ] - X B ) ) , m e a n = c ( 0 , 0 ) , s i g m a = n e g s i g m a ) - p m v n o r m ( l o w e r = c ( - I n f , - I n f ) ,   u p p e r = c ( Z G , ( t a u [ 1 ] - X B ) ) , m e a n = c ( 0 , 0 ) , s i g m a = n e g s i g m a ) ) }  
 	 	 e l s e   i f ( i d e n t i c a l ( Y [ i ] , 2 ) ) { l l i k [ i ] < - l o g ( p m v n o r m ( l o w e r = c ( - I n f , - I n f ) ,   u p p e r = c ( Z G , ( t a u [ 3 ] - X B ) ) , m e a n = c ( 0 , 0 ) , s i g m a = n e g s i g m a ) - p m v n o r m ( l o w e r = c ( - I n f , - I n f ) ,   u p p e r = c ( Z G , ( t a u [ 2 ] - X B ) ) , m e a n = c ( 0 , 0 ) , s i g m a = n e g s i g m a ) ) }  
 	 	 e l s e   i f ( i d e n t i c a l ( Y [ i ] , 3 ) ) { l l i k [ i ] < - l o g ( p m v n o r m ( l o w e r = c ( - I n f , - I n f ) ,   u p p e r = c ( Z G , ( X B - t a u [ 3 ] ) ) , m e a n = c ( 0 , 0 ) , s i g m a = s i g m a ) ) }            
 	 	 }  
 	 l l i k < - - 1 * s u m ( l l i k ) 	 	 	 	 	              
 	 r e t u r n ( l l i k )  
 	 } 	 	 	 	 	 	 	 	              
  
  
  
 # s e t   s t a r t i n g   p a r a m e t e r s  
 e s t < - r b i n d ( . 1 , . 1 , . 1 , . 1 , . 1 , . 1 , . 1 , . 1 , . 1 )  
  
 # s e t   d a t a ,   Y ,   X ,   a n d   Z  
 d a t a < - d a t a s e t  
 Y < - d a t a s e t $ y  
 X < - c b i n d ( d a t a s e t $ x 1 , d a t a s e t $ x 2 )  
 Z < - c b i n d ( 1 ,   d a t a s e t $ x 1 ,   d a t a s e t $ x 2 )  
  
 # o p t i m i z e  
 o u t p u t . Z I O P C < - n l m ( f = Z I O P C ,     p = e s t ,   Y = Y , X = X , Z = Z   , i t e r l i m = 5 0 0 ,   d a t a = d a t a s e t ,   h e s s i a n = T R U E )  
 p r i n t ( o u t p u t . Z I O P C )  
 v c v < - s o l v e ( o u t p u t . Z I O P C $ h e s s i a n )  
 v c v    
  
  
 